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Abstract 

The purpose of this study is to propose a high-accuracy and fast numerical method 
for the Cauchy problem of the Laplace equation. Our problem is directly dis- 
cretized by the method of fundamental solutions (MFS). The Tikhonov regular- 
ization method stabilizes a numerical solution of the problem for given Cauchy data 
with high noises. The accuracy of the numerical solution depends on a regulariza- 
tion parameter of the Tikhonov regularization technique and some parameters of 
MFS. The L-curve determines a suitable regularization parameter for obtaining an 
accurate solution. Numerical experiments show that such a suitable regularization 
parameter coincides with the optimal one. Moreover, a better choice of the parame- 
ters of MFS is numerically observed. It is noteworthy that a problem whose solution 
has singular points can successfully be solved. It is concluded that the numerical 
method proposed in this paper is effective for a problem with an irregular domain, 
singular points, and the Cauchy data with high noises. 
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1 Introduction 



Many kinds of inverse problems have recently been studied in science and en- 
gineering. The Cauchy problem of an elliptic partial differential equation is a 
well known inverse problem. The Cauchy problem of the Laplace equation is 
an important problem which can be applied to inverse problem of electrocar- 
diography [12]. Onishi et al. [9] proposed an iterative method for solving the 
Cauchy problem of the Laplace equation. This method reduces the original 
inverse problem to an iterative process which alternatively solves two direct 
problems. This method, called the adjoint method in the papers [7], [10], can 
solve various inverse problems by applying many kinds of numerical methods 
for solving partial differential equations, such as the finite difference method 
(FDM), the finite element method (FEM), and the boundary element method 
(BEM). The convergence of this method for the Cauchy problem of the Laplace 
equation has been obtained [6]. 

The method of fundamental solutions (MFS) is effective for easily and rapidly 
solving the elliptic well-posed direct problems in complicated domains. Mathon 
and Johnston [13] first showed numerical results obtained by the MFS. The 
papers [11], [14] discuss some mathematical theories on the MFS. Both of 
the BEM and the MFS are well known boundary methods, which discretize 
original problems based on the fundamental solutions. The MFS does not 
require any treatments for the singularity of the fundamental solution, while 
the BEM requires singular integrals. The MFS is a true meshless method, and 
can easily be extended to higher dimensional cases. 

Wei et al. [4] applied the MFS to the Cauchy problems of elliptic equations. 
This method uses the source points distributed outside the domain. The ac- 
curacy of numerical solutions depends on the location of the source points. 
They numerically showed the relation between the accuracy and the radius of 
a circle where the source points are distributed. But, the relation between the 
accuracy and the number of source points has not clearly been given, yet. 

Many researchers have solved the Cauchy problem by various methods. How- 
ever, to our knowledge, the conventional methods cannot solve a problem 
whose solution has singular points (see [8] for example). 

In this paper, we use the MFS to directly discretize the Cauchy problem of 
the Laplace equation. This problem is an ill-posed problem, where the solution 
has no continuous dependence on the boundary data. Namely, a small noise 
contained in the given Cauchy data has a possibility to affect sensitively on 
the accuracy of the solution. The problem is discretized directly by the MFS 
and an ill-conditioned matrix equation is obtained. A numerical solution of the 
ill-conditioned equation is unstable. The singular value decomposition (SVD) 
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can give an acceptable solution to such an ill-conditioned matrix equation. The 
SVD was successfully applied to the MFS for solving a direct problem [15]. 
Even though we apply the SVD, we still cannot obtain an acceptable solution 
for the case of the noisy Cauchy data. We use the Tikhonov regularization 
to obtain a stable regularized solution of the ill-conditioned equation. The 
regularized solution depends on a regularization parameter. Then, we need to 
determine a suitable regularization parameter to obtain a better regularized 
solution. Hansen [1] suggested the L-curve as a method for finding the suitable 
regularization parameter. It is known that the suitable parameter is the one 
corresponding to a regularized solution near the "corner" of the L-curve. We 
can find the corner of the L-curve as a point with the maximum curvature [5]. 

Under the assumption of uniform distribution of the source and the colloca- 
tion points, we will numerically indicate that a suitable regularized solution 
obtained by the L-curve is optimal in the sense that the error is minimized. We 
will respectively show the accuracy and the optimal regularization parameter 
against a noise level. We will also mention influence of the total numbers of the 
source and the collocation points on accuracy. We will show that our method 
is effective for a problem whose solution has singular points. It is noteworthy 
that such kind of problems can also successfully be solved. 

Section 2 introduces the Cauchy problem. In Section 3, the MFS discretizes 
the problem. In Section 4, the singular value decomposition, the Tikhonov reg- 
ularization and the L-curve are used to obtain a suitable regularized solution. 
In Section 5, numerical experiments confirm that the suitable regularization 
parameter by the L-curve coincides with the optimal one that minimizes the 
error between the regularized solution and the exact one. The error and the 
optimal regularization parameter against the noise level of the Cauchy data 
are respectively shown. Then, our interest is how to choose the following three 
parameters used in MFS: the numbers of collocation points, the number of 
source points, and the radius of a circle where source points are distributed. 
A better choice of the parameters is also observed. A problem with an irregu- 
lar domain and a problem whose solution has singular points are successfully 
solved, respectively. Section 6 concludes the paper. 



2 Problem Setting 



We consider the Laplace equation — Au = in a two-dimensional bounded 
domain fl enclosed by the boundary T. We prescribe Dirichlet and Neumann 
boundary conditions simultaneously on a part of the boundary T, denoted by 
r 1; as follows: 

, du 
u = f, tt = 9 on ri, 
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where / and g denote given continuous functions defined on r\, and n the 
unit outward normal to r\. Then, we need to find the boundary value u on 
the rest of the boundary T 2 '■— T \ F 1 or the potential u in the domain Q. 
This problem is called the Cauchy problem of the Laplace equation, and the 
boundary data are called the Cauchy data. 

Our Cauchy problem is described as follows: 



Problem 1 For the given Cauchy data f,g e C(Ti), find u e C(T 2 ) or 
u e C 2 {Vt) n C^U) such that 



-Au = in n, (1) 
« = /» 9^ = 9 on r\. (2) 



The Cauchy problem is a well known ill-posed problem. We can show the 
instability of the solution to the Cauchy problem of the Laplace equation as 
follows: For example, in the case where 



n=(0,l) 2 = {(x,y):0<x<l, < y < 1}, 
T = [0, 1] x {0} = {(x, 0) : < x < 1}, 

f(x,0) = -\sm(nx), g(x,0)=0 (k > 0), 



n k 



the solution is given by 



u(x, y) = —r sin(nx) cosh(ra/). 



n k 



Here, we can see that 



sup |/(a?)| — > 0, sup |m(£c)| — > oo, 



from which we know that the solution u of the Cauchy problem does not 
depend continuously on the Cauchy data / and g. 
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3 Discretization by the Method of Fundamental Solutions 




&N-1 



Fig. 1. Problem setting and distributions of collocation and source points 



The fundamental solution of the Laplace equation in two dimensions is denned 

as 

ip*(r) := lnr 



for r = | a; | = ^Jx 1 + y' 2 , which is a solution to 

-Ay?* = 5(x). 



We distribute the collocation points {xi\^L l C I\ on the boundary where the 
Cauchy data are prescribed, and the source points {^ J }^ 1 C Q along a circle 
outside the domain (Figure 1). We approximate u by ujy: 



N 

u{x) ~ Un(x) := Wjipj(x), 

3=1 



(3) 



where the basis function is defined as 



(Pj(x) :=(p*{\x-SA) 



(4) 



and {vjj}j =1 are expansion coefficients to be determined below. Since the ba- 
sis functions (4) have no singular points in fl, the approximate function un 
satisfies the Laplace equation (1). Substituting (3) into (2) and assuming that 
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(2) is satisfied at the collocation points, we have 



( N 



•■ .7 = 1 



or in the matrix form: 



Aw = b, 



where the matrix A = (a^) G R 2MxN an d the vectors w 
b = (bi)e R 2M are defined by 



(5) 



dip j 



i = l,2,...,M 



dn 



x, A/ ). i = M + 1, M + 2, . . . , 2M 



j = l,2,. ..,N, 



hi := < 



/(xO, t = l,2,...,M 
g(xi- M ), i = M + 1, M + 2, . . . , 2M 



4 Regularized Solution 

4-1 Singular value decomposition 

In general, (5) has the case where the exact solution w does not exist in the 
conventional sense. As an exact solution, we consider the least square and 
least norm solution w defined by 

||u>o|| = min llioll, W := {w : \\Aw — b\\ = min \\Aw — 611}. (6) 

w<=W w£R N 

In this paper, we refer the solution w as the exact solution to (5). 

For the matrix A e R 2MxN , the singular value decomposition (SVD) can be 
written as follows: 

A = UnV T = ^2(T i u i vf, 

i=l 
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where 

TT („. „. \ r- r>2Mx2M i; /„. „. „. \ ,- r>NxN 

u — (iii« 2 • • • u 2 m) e« , v = {vi v 2 ■ ■ ■ v N ) e k 



i 



Oi (i = 3, 



E = ( % ) G R 2MxN , s t 

" 



UU T = I 2M e R 2M * 2M , VV T = I N e R NxN , 
2M > N, o x > ■ ■ ■ o r > 0, o r+1 = • • • = o N = 



with the identity matrices J2M and 7jv. The non-negative values {<Ti}^Li are 
called the singular values of the matrix A. Using the SVD, we can express the 
solution to (6) as 

i=i a * 

In a real problem, the Cauchy data / and 5 contain some noises. We consider 
the following equation instead of (5): 

Aw = b s , b 5 = b + Ab (7) 

with the noise vector Ab e -R 2M . Then, the solution w s to 

\\w s \\ = min \\w\\, W 5 :={w:\\Aw-b 5 \\= min \\Aw-b 5 \\} (8) 
mew* w<=r n 

is quite different from the exact solution w to (6) since the solution is dis- 
continuous for the Cauchy data. We need to find a good approximation to 
w . 



4-2 Tikhonov regularization 



In order to obtain a good approximate solution to (8), we consider minimizing 
the following functional with a regularization parameter a > according to 
the Tikhonov regularization: 

J s a (w) := \\Aw-b s \\ 2 + a 2 \\w\\ 2 . (9) 

It is easy to see that the functional is strictly convex for any a > 0. Hence, 
J s a has a unique minimum point w s a called the regularized solution: 

■£(«£)= min 

We know that is the solution to 

(A T A + a 2 I N )w s a = A T b s . (10) 
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The equation (10) is uniquely solvable since the matrix (A T A + a 2 lN) is sym- 
metric positive definite. The SVD of w 5 a can be expressed in the form: 

s * (u u b s ) 

w « = Z^li—Z V i 

i=l a i 

with the filter factor % := of /(of + a 2 ). Then, substituting w 5 a into (3), we 
find the approximate potential Un in VL U T 2 . 

The error between the regularized solution for the noisy data and the exact 
solution is decomposed into 

wi- Wo = K - O + « - « ) = E 7.^^ + E(7* - 1)^^, 

1=1 a * 1=1 a * 

(11) 

The first term is the perturbation error due to the relative noise Ab and the 
second term is the regularization error caused by regularization of the exact 
b. When < a < 1, we see that 7$ ~ 1 for most of i, and the error w s a — w 
is dominated by the perturbation error. On the other hand, when a 1, we 
see that 7^ <^ 1 and the error w s a — w is dominated by the regularization 
error. In the next subsection, we will consider a useful method for finding a 
suitable regularization parameter to minimize both of the perturbation and 
the regularization errors. 



4-3 L- curve 

To find a suitable regularization parameter, Hansen [1] suggests the L-curve, 
which is defined as the continuous curve consisting of all the point (|| — 
b*||,||u£||) for a >0: 

£:={(\\Awi-b s \\,\\w s J):a>0}. 

For fixed a > 0, we get w 5 a and then can calculate the residual norm H^w^ — 
b s \\ and the solution norm Thus, the L-curve can be plotted as a set 

of all the points of the residual norms as abscissa and the solution norms as 
ordinate for all a > 0. 

The L-curve is plotted in double logarithm, and displays the compromise be- 
tween minimization of the perturbation error and the regularization error in 
(11). A suitable regularization parameter is given by the one corresponding to 
a regularized solution near the "corner" of the L-curve. The "corner" can be 
regarded as the point where the curvature of the L-curve becomes maximum 
[3], [5]. 
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In the case when 2M = N and the Cauchy data /, g have no noises, if A is 
non-singular, we can directly solve (5) to obtain a solution with high accu- 
racy. However, we cannot guarantee that (5) is always solvable. Even if there 
exists the inverse matrix, the solution to (5) for the noisy Cauchy data differs 
from the exact solution. On the other hand, the regularized solution by the 
Tikhonov regularization is always uniquely determined for a > 0. In the next 
section, our numerical experiments will show that the suitable regularization 
parameter given by L-curve coincides with the optimal one a opt defined by 

\\w -w 5 a II = min \\w - w s J\. 

II u a op t H a>Q H an 



5 Numerical Experiments 



5.1 Circular domain 



We first consider a harmonic function u(x,y) = e x cosy — e y smx in a unit 
disk VL := {(x, y) : x 2 + y 2 < 1}. According to the exact potential u, the exact 
Cauchy data are given by / = u and g = du/dn on the fourth part of the 
whole boundary T, which is defined by 

I\ :={{x,y) :x 2 + y 2 = l, x > 0, y > 0}. 

We now assume that the exact potential u is unknown, and identify a boundary 
value on the rest of the boundary T 2 :— dQ \ Ti from the noisy Cauchy data 
f s = (1 + s)f and g 5 = (1 + e)g, where e = e(x, y) is a uniform random noise 
such that —5 < e(x,y) < 5 with the relative noise level of 1005%. 



We distribute uniformly the collocation points {a^}^ C T 1 and the source 



points {£j}f =1 C fi c as follows: 



27[(i — 1) 7T 

Xi = (cos 0;, sin 0;), 9i = — — — + jtj, i = l,2,...,M, 

- 2ir(i -1) it { ' 

^ = (Rcos9 3 ,Rsme 3 ), 6j = ^ ' + -, j = 1, 2, . . . , N, 

where R > 1 is the radius of the circle where the source points are distributed. 
We adopt the MATLAB code for solving discrete ill-posed problems based on 
SVD, made by Hansen [2], [3], to our numerical computations. Due to the 
maximum principle, it is sufficient to confirm the boundary error between 
the identified potential un and the exact one u rather than the domain error 
in our numerical experiments. We define the maximum relative error on the 
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boundary by 



e := 



\Un — u\ 



u 



where the maximum norm on the boundary denotes 



|w||oo = sup |m(je)|, \/u £ C(T). 




-4-3-2-1 1 2 3 



Fig. 2. Distributions of collocation and source points ((R,M,N) = (3.2,600,28)) 
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Fig. 3. L-curve (The corner is located at the point for a = 2.1195 x 10 ) 
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Fig. 4. The regularized solutions versus the exact solution (5% noise level) 



Table 1 

The maximum relative errors for each a 



a 


2.12 x 1(T 4 


2.12 x 1(T 3 


2.12 x lO" 2 


error 


0.2743 


0.0305 


0.3382 



In the first experiment, the relative noise level of the Cauchy data is assumed to 
be 5% (5 = 0.05). We set the parameters (R, M, N) = (3.2,600,28). Figure 2 
shows the distributions of the collocation and the source points. As we can 
see in Figure 3, the corner of the L-curve is located at the point (||y4iy^ — 
b s \\, || 10* ||) with the regularization parameter a = 2.1195 x 10~ 3 . Figure 4 
shows the regularized solutions on the boundary for a = 2.1195x 10 -4 , 2.1195x 
10 -3 , 2.1195 x 10~ 2 , 0. We can see that the solution is quite unstable if a = 0, 
that is, if the regularization is unapplied. Comparing the other solutions for 
a = 2.1195 x 10" 4 ,a = 2.1195 x 10" 3 , 2.1195 x 10~ 2 , we can confirm that 
a = 2.1195 x 10~ 3 is a suitable regularization parameter to obtain a better 
approximate solution (Table 1). 
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a 

Fig. 5. Maximum relative error against regularization parameter 

From Figure 5, we can see that the maximum relative error reaches a minimum 
at a = io~ 2 - 673 2.12 x 10~ 3 , which coincides with the suitable regularization 
parameter obtained by the L-curve. Hence, we know that the optimal regu- 
larization parameter can be given as the one corresponding to a regularized 
solution at the corner of the L-curve. 




-14 -12 -10 -8 -6 -4 -2 -14 -12 -10 -8 -6 -4 -2 

logio * lo gio s 

(a) (b) 

Fig. 6. (a) Maximum relative error and (b) the optimal regularization parameter 
against relative noise level 

Figure 6 (a) shows the maximum relative error e for the optimal regularized 
potential against the relative noise level S. The regression line in the interval 
[—9,0] is expressed by log 10 e = 0.37951 log 10 5 — 0.22672. For the optimal 
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regularized potential un, we have e = 0(S 0,38 ) for 5 > 1CT 9 . Figure 6 (b) 
indicates the optimal regularization parameter a opt against the relative noise 
level 6 and the regression line in the interval [—9,0] given by log 10 a opt = 
1.0186 log 10 5 — 1.2434. From this numerical result, we can obtain the relation 
« op t = 0(5) for 10- 9 < 5 < 1. 

After setting the parameters (R,M,N), we can obtain a suitable regularized 
solution based on the Tikhonov regularization and the L-curve. Now, our 
problem is how to choose suitable parameters (R, M, N). 

Figure 7 shows the maximum relative error against the number of colloca- 
tion points. We know from this result that we need to take sufficiently many 
collocation points to obtain accurate solutions. 




o ' 1 1 1 1 1 1 

500 1000 1500 2000 2500 3000 

M 



Fig. 7. Maximum relative error against number of collocation points 




N N 

(a) (b) 

Fig. 8. Contour line of the maximum relative error against (N,R) (M = 600); (a) 
(N,R) e [10,60] x [2,12], (b) (N,R) & [24,30] x [3,3.4] 
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Figure 8 shows the contour line of the maximum relative error e against (N, R) 
for the fixed number of collocation points M = 600. Through this result, we 
know that the maximum relative error is roughly independent of the number 
of source points N for the fixed radius R of the circle where source points 
are distributed, and becomes large for large R. As a result, we know that the 
parameters R ~ 3.2 and N > 25 will yield a better regularized solution. 



5.2 Irregular domain 



As the next example, we assume the exact solution as same as the one in the 
previous example in an irregular domain enclosed by the boundary 

x{6) = r(9) cos 6, y(6) = r(6) sin < 6 < 2vr 

with the Cassini oval 



r(0) = r(6; a, b) = a\j cos 26 + y/ {b/a) 4 - sin 2 26, < 6 < 2tt (13) 

in the polar coordinates with a = 1, b = 1.01. It is easy to show that the unit 
outward normal to the boundary is expressed as 



n{6) 



r'(6) sin 6 + r(6) cos 6 -r'{9) cos 6 + r(6) sin 6 
^(r'(6)) 2 + r(6) 2 ' yj(r'(6)) 2 + r(9) 2 



We distribute collocation and source points along the boundary and the circle 
(R cos 6 , R sin 6) uniformly as similar as (12). Figure 9 shows the domain, the 
unit outward normal, the collocation and the source points for example. 



i - + 




Fig. 9. Domain and distribution of points ((R, M, N) = (2, 40, 10)) 
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The relative noise level of the Cauchy data is assumed to be 10% (5 = 0.1). Let 
(R, M, N) = (2, 5200, 30). The optimal regularization parameter can be found 
as a = 2.3493 x 10~ 2 by using the L-curve. Figure 10 shows the regularized 
solution on the boundary with respect to the optimal regularization parameter. 

3.5 

3 
2.5 

2 

u 

1.5 



0.5 

jt/2 tt 3tt/2 2tt 

6 

Fig. 10. The regularized solution ujy versus the exact u (10% noise level; 
(R,M,N) = (2,5200,30)) 

From the result, it is concluded that even if the boundary of the domain is 
complicated and the noise level of the Cauchy data is higher, the regularized 
solution is in very good agreement with the exact one. 

5.3 Problems with singular points 

We consider two problems whose solutions have singular points. 
We first assume the exact solution 

u(x, y) = log y/(x- 0.2) 2 + y 2 - log y/(x + 0.2) 2 + y 2 
in the annulus domain 

n = {(x,y) : 0.5 2 < x 2 + y 2 < 1} 
with the outer and the inner boundaries 

Tout = {(x, y):x 2 + y 2 = 1}, r in = {(x, y) : x 2 + y 2 = 0.5 2 }. 
The exact solution u has two singular points at (x,y) = (—0.2,0), (0.2,0). 
We now assume that the exact potential u is unknown. From the Cauchy data 



regularized 
exact 
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given on the fourth part of the outer boundary r out , defined by 



Ti := {(x,y) : x 2 + y 2 = 1, x > 0, y > 0}, 

we identify a boundary value on the rest of the boundary (T out \ Ti) U r in . 

We distribute uniformly the collocation points {xi}fi ± C Fx- The source points 
C Q are uniformly distributed along two circles whose centers are 
and radii are R out and Ri n , respectively. 

In the first case, the exact Cauchy data is assumed to be given. Let (-R ou t, -Rin, M, 
(3.2,0.4,600,60) (Figure 11). Figure 12 shows the identified solution on the 
outer and the inner boundaries. We can see that the identified solution is in 
very good agreement with the exact one in spite of the solution with singular 
points. 




-4-3-2-1 1 2 3 4 



Fig. 11. Domain and distribution of points ((-Rout, -Rin, M, N) = (3.2,0.4,600,60)) 
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vr/2 7T 3vr/2 2tt 
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Fig. 12. The regularized solution un versus the exact u (no noise) 
In the second case, the relative noise level of the Cauchy data is assumed to 
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identified (outer) 




be 5% (5 = 0.05). Let (R out , R in , M, N) = (3.2,0.2,5200,30) (Figure 13). The 
optimal regularization parameter can be found as a = 6.3932 x 10 -3 by using 
the L-curve. Figure 14 shows the regularized solution on the outer and the 
inner boundaries with respect to the optimal regularization parameter. From 
this result, we know that the regularized solution is acceptable. 



* 




Fig. 13. Domain and distribution of points ((-Rout, -Kin, M, N) = (3.2,0.2,5200,30)) 
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Fig. 14. The regularized solution un versus the exact u (5% noise level) 

As another example, we assume that the exact solution is given by 

u(x,y) = X 

x 2 + y l 

in the same domain as above. The Cauchy data with 5% noise level are pre- 
scribed on the same part of the boundary as above. Let (R ont , Ri n , M, N) = 
(3.2, 0.05, 5200, 30). Figure 15 shows the regularized solution on the outer and 
the inner boundaries with respect to the optimal regularization parameter 
a = 1.5088 x 10 -3 . The accuracy of the regularized solution is quite good. 
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vr/2 tt 3tt/2 



Fig. 15. The regularized solution un versus the exact u (5% noise level) 



6 Conclusions 



We consider using MFS as the numerical method for the Cauchy problem of 
the Laplace equation. Since MFS is a messless method, we can easily treat 
a complicated boundary. This paper proposes a direct method instead of an 
iterative one. The Tikhonov regularization can find a stable solution. The L- 
curve automatically gives a suitable regularization parameter, which coincides 
with the optimal one shown in the numerical experiments. Hence, after setting 
the parameters (R, M, N) the optimal regularized solution can be obtained 
quickly and automatically. Moreover, our numerical method can successfully 
solve even a problem whose solution has singular points. 

The following is the guideline for choosing better parameters (R, M, N): The 
collocation points should be distributed as many as possible compared with 
the source points. There is no value increasing the number of source points 
N. It is enough for iV « 30 to obtain a better solution. The radius of the 
circle where source points are distributed should be small like R ~ 3, since 
the stability is more important than the accuracy in this inverse problem. 

In conclusion, the numerical method proposed in this paper is applicable for 
solving a problem in a complicated domain with the Cauchy data that contains 
large noises even with a noise level of 10%. This method is also effective for 
solving even a problem with singular points. 
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